با توجه به داده بانک جهانی به سوالات زیر پاسخ دهید. برای استفاده از داده از سه فایل زیر استفاده نمایید. داده نام کشورها: WDICountry داده نام سری های زمانی: WDISeries داده کل: WDIData در صورتی که داده را در اختیار ندارید می توانید از بسته WDI استفاده نموده و داده های مورد نظر را استخراج نمایید.
بارگذاری داده ها و کتابخانه ها:
library(readr)
library(dplyr)
library(ggplot2)
library(highcharter)
library(stringr)
wdi_data = read_csv('data/WDIData.csv')
۱. ده کشور فقیر دنیا را بیابید. نمودار درآمد روزانه آنها را رسم کنید. چند درصد از اعضای این کشورها زیر خط فقر هستند؟ متوسط عمر در این کشورها چقدر است؟
برای حل این سوال ابتدا ستون های مورد نظر را برای سال ۲۰۱۶ پیدا می کنیم. فقط از آنجایی که خط فقر NA زیادی داشت از سال ۲۰۱۴ آن استفاده می کنیم. سپس بر اساس فقر مرتب کرده و ده کشور فقیر را انتخاب می کنیم. سپس خط فقر و امید به زندگی را رسم می کنیم. در نهایت نمودار را رسم می کنیم.
poverty <- wdi_data %>% filter(`Indicator Code`== 'NY.GDP.PCAP.PP.CD') %>%
select(country = `Country Name`, poverty = `2016`)
income <- wdi_data %>% filter(`Indicator Code`== 'NY.ADJ.NNTY.CD') %>%
select(country = `Country Name`, income = `2016`)
poverty_line <- wdi_data %>% filter(`Indicator Code`== 'SI.POV.NAHC') %>%
select(country = `Country Name`, poverty_line = `2014`)
life_expectancy <- wdi_data %>% filter(`Indicator Code`== 'SP.DYN.LE00.IN') %>%
select(country = `Country Name`, life_expectancy = `2016`)
population <- wdi_data %>% filter(`Indicator Code`== 'SP.POP.TOTL') %>%
select(country = `Country Name`, population = `2016`)
poverty <- poverty %>% inner_join(income, by = c("country")) %>%
inner_join(population, by = c("country")) %>%
inner_join(poverty_line, by = c("country")) %>%
inner_join(life_expectancy, by = c("country")) %>%
arrange(poverty) %>%
mutate(daily_income = income/(365*population)) %>%
slice(1:10)
knitr::kable(poverty)| country | poverty | income | population | poverty_line | life_expectancy | daily_income |
|---|---|---|---|---|---|---|
| Central African Republic | 698.7067 | 1374375431 | 4594621 | NA | 52.171 | 0.8195262 |
| Burundi | 777.7529 | 2122720795 | 10524117 | 64.9 | 57.481 | 0.5526044 |
| Congo, Dem. Rep. | 801.6301 | 26138587178 | 78736153 | NA | 59.621 | 0.9095259 |
| Liberia | 812.6739 | 1112973037 | 4613823 | 54.1 | 62.505 | 0.6608925 |
| Niger | 986.2070 | 6041624602 | 20672987 | 44.5 | 60.058 | 0.8006775 |
| Malawi | 1168.8256 | 3635417843 | 18091575 | NA | 63.223 | 0.5505352 |
| Mozambique | 1216.7928 | 8566275076 | 28829476 | 46.1 | 58.311 | 0.8140712 |
| Sierra Leone | 1476.2137 | 2676603250 | 7396190 | NA | 51.835 | 0.9914780 |
| Togo | 1490.5362 | 3020251309 | 7606374 | NA | 60.232 | 1.0878588 |
| Madagascar | 1506.2383 | 8910175309 | 24894551 | NA | 65.932 | 0.9805937 |
poverty %>% arrange(daily_income) %>%
hchart(type = "column", hcaes(x = country, y = daily_income)) %>%
hc_yAxis(title = list(text = "Daily Income")) %>%
hc_xAxis(title = list(text = "Country")) %>%
hc_title(text = "Daily Income in Poor Countries", style = list(fontWeight = "bold")) %>%
hc_add_theme(hc_theme_ffx())۲. تراژدی روآندا: بر اساس داده های امید به زندگی ابتدا نمودار سالانه نمودار جعبه ایی امید به زندگی کشورها را رسم نمایید(در یک نمودار!). سپس سری زمانی امید به زندگی روآندا را به آن اضافه کنید. چه می بینید؟ چند میلیون نفر کشته شدند؟
برای حل این سوال، امید به زندگی را در تمامی سال ها بدست می آوریم. سپس نمودار را از ستونی به ردیفی به کمک melt تبدیل می کنیم. در نهایت نمودار جعبه را می کشیم.
library(reshape2)
life_expectancy <- wdi_data %>% filter(`Indicator Code`== 'SP.DYN.LE00.IN') %>%
select(country = `Country Name`, matches('\\d{4}')) %>%
melt(id.vars=c('country'))
rwd_life_expectancy <- life_expectancy %>% filter(country == 'Rwanda')p = ggplot(data = life_expectancy, mapping = aes(x = variable, y = value)) + geom_boxplot() +
geom_line(data = rwd_life_expectancy, mapping = aes(x = variable, y = value, group = 1, color = 'Rwanda')) +
xlab("Year") + ylab("Life Expectancy") + ggtitle("Average Life Expectancy Worldwide") +
theme(axis.text.x = element_text(angle = 90, hjust = 1))
p
مشاهده می کنیم که در این سال ها امید به زندگی در روآندا با اختلاف از سایر کشورها کمتر است که علت آن نسل کشی توسط دولت بوده است. همچنین در زیر می بینیم که یک میلیون نفر در این دوران کشته شده اند.
rwd_death <- wdi_data %>% filter(`Indicator Code`== 'SP.DYN.CDRT.IN') %>%
select(country = `Country Name`, matches('\\d{4}')) %>% filter(country == 'Rwanda') %>%
melt(id.vars=c('country')) %>% filter(!is.na(value)) %>%
summarise(tot_death = 1000*sum(value))
cat("Total Death in Rwanda is", rwd_death[1,], "people.")Total Death in Rwanda is 1053185 people.۳. نمودار امید به زندگی و هزینه های بهداشتی را رسم کنید. چه نتیجه ایی می گیرید؟
برای حل این سوال همانند بالا داده را در باره ی هزینه های بهداشتی نیز بدست می آوریم، منتهی برای آنکه رنج آن با امید به زندگی قابل مقایسه باشد داده ها را به ۱۰۰ تقسیم می کنیم و دو نمودار را کنار یکدیگر رسم می کنیم.
health <- wdi_data %>% filter(`Indicator Code`== 'SH.XPD.CHEX.PC.CD') %>%
select(country = `Country Name`, matches('\\d{4}')) %>%
melt(id.vars=c('country')) %>%
mutate(value = value/100, type = 'health') %>%
filter(!is.na(value))
life_expectancy <- life_expectancy %>% filter(!is.na(value)) %>% mutate(type = 'life expectancy')
health_expect <- health %>% rbind(life_expectancy)
p = ggplot(data = health_expect, mapping = aes(x = variable, y = value, fill = type, color = type)) + geom_boxplot() +
xlab("Year") + ylab("Life Expectancy or Health Cost") + ggtitle("Average Life Expectancy and Health Cost Worldwide") +
theme(axis.text.x = element_text(angle = 90, hjust = 1)) 
همانطور که انتظار می رود، هر چه امید به زندگی به طور میانگین بالا می رود، میزان هزینه های بهداشتی و اهمیت افراد به بهداشت خود نیز افزایش پیدا می کند.
۴. آیا قدرت خرید خانواده های ایرانی در ۵۰ سال اخیر افزایش یافته است؟ برای این کار از داده های اقتصادی خانوار استفاده کنید.
برای این سوال همانند سوال بالا ردیف مناسب با ستون های مناسب را به ازای تمامی سال ها برای ایران بدست می آوریم و نمودار را می کشیم.
iran_buy_cap <- wdi_data %>% filter(`Indicator Code`== 'NY.GDP.PCAP.PP.CD') %>%
select(country = `Country Name`, matches('\\d{4}')) %>%
filter(country == 'Iran, Islamic Rep.') %>%
melt(id.vars=c('country')) %>%
filter(!is.na(value))
iran_buy_cap %>% arrange(variable) %>%
hchart(type = "line", hcaes(x = variable, y = value)) %>%
hc_yAxis(title = list(text = "Purchasing Capability")) %>%
hc_xAxis(title = list(text = "Year")) %>%
hc_title(text = "Purchasing Capability in Iran", style = list(fontWeight = "bold")) %>%
hc_add_theme(hc_theme_economist())همانطور که می بینیم قدرت خرید در ایران زیاد شده است، اما در پایان سال ۲۰۱۲ که با پایان دولت قبلی همزمان است، میزان قدرت خرید به علت اختلاس ها و عدم کنترل وضعیت اقتصادی افت کرده است.
۵. رشد اقتصادی ایران را با کشورهای دیگر در طول ۲۰ سال گذشته بر حسب بیست شاخص های اقتصادی مهم مانند تولید ناخالص ملی، تورم و … ارزیابی کنید! (برای هر شاخص از تصویرسازی استفاده کنید.)
برای حل این سوال ۲۰ معیاری که مقادیر NA کمتری دارد را می یابیم و در حل مرحله داده ی جهانی را با ایران مقایسه می کنیم. نام شاخص در عنوان خطوط نمودار آمده است. البته مشاهده می کنیم که دو نمودار تقریبا داده ای ندارند.
world_ec <- wdi_data %>% filter(`Indicator Code`== 'NE.EXP.GNFS.CD') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]')) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'NY.GDP.PCAP.CD') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]'))) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'NY.GNP.MKTP.CD') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]'))) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'NE.GDI.FTOT.KD.ZG') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]'))) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'BX.KLT.DINV.CD.WD') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]'))) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'BM.GSR.GNFS.CD') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]'))) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'BN.GSR.FCTY.CD') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]'))) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'FP.CPI.TOTL.ZG') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]'))) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'FM.LBL.BMNY.GD.ZS') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]'))) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'NV.AGR.TOTL.ZS') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]'))) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'NE.CON.TETC.ZS') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]'))) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'NE.CON.GOVT.ZS') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]'))) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'NV.IND.MANF.ZS') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]'))) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'MS.MIL.XPND.GD.ZS') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]'))) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'TG.VAL.TOTL.GD.ZS') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]'))) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'NY.GDP.TOTL.RT.ZS') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]'))) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'NE.TRD.GNFS.ZS') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]'))) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'FS.AST.DOMS.GD.ZS') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]'))) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'NY.GNP.MKTP.PP.CD') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]'))) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'SH.UHC.NOP1.TO') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]')))
world_ec_m <- world_ec %>%
melt(id.vars=c('country', 'index')) %>%
filter(!is.na(value)) %>%
group_by(index)
plots <- world_ec_m %>%
do(
plots = ggplot(.,aes(x = variable, y = value, color = index)) + geom_boxplot() +
xlab("Year") + ggtitle("Worldwide Economic vs. Iran") +
theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
geom_line(data = subset(.,country == 'Iran, Islamic Rep.'),aes(x = variable, y = value, group = 1, color = 'Iran, Islamic Rep.'))
)
plots$plots[[1]]
[[2]]
[[3]]
[[4]]
[[5]]
[[6]]
[[7]]
[[8]]
[[9]]
[[10]]
[[11]]
[[12]]
[[13]]
[[14]]
[[15]]
[[16]]
[[17]]
[[18]]
[[19]]
[[20]]
همانطور که در نمودارها می بینیم، ایران تقریبا بر روی میانگین کشورها قرار دارد، اما مقادیری همچون هزینه ی مصرف نهایی از کشورهای دیگر بسیار کمتر است و سرمایه گذاری در خارج آن صفر است. هم چنین میزان تجارت از میانگین کمتر است و میزان اجاره ی منابع طبیعی بسیار بیشتر از سایر کشورها است. در کل ایران اقتصادی متوسط رو به پایین دارد.
۶. در قسمت قبل با استفاده از روش خوشه بندی k-means داده ها را به سه دسته تقسیم کنید. ایران در کدام دسته می گنجد؟ (پیش از خوشه بندی طبیعتا داده را باید پاکسازی و استاندارد سازی نمایید.)
برای پاکسازی داده ها مقادیر ستون های NA را برابر با صفر قرار می دهیم، سپس برای هر کشور، میانگین مقادیر ستون ها را محاسبه می کنیم. در نهایت به کمک kmeans داده های سال های مختلف برای کشورهای مختلف را به سه دسته تقسیم می کنیم.
# cleaning data
world_ec_cluster <- world_ec
world_ec_cluster[is.na(world_ec)] = 0
patt = "Africa|America|Asia|dividend|Euro|income|members|Middle|only|Sub|total|World"
world_ec_cluster <- world_ec_cluster %>% group_by(country) %>%
filter(!grepl(patt,country)) %>%
summarise_at(c(2:22), mean)
set.seed(1234)
clusters = kmeans(world_ec_cluster[,2:22], 3)
world_ec_cluster$cluster_no = as.integer(clusters$cluster)
iran_cluster = world_ec_cluster %>%
filter(country == 'Iran, Islamic Rep.')
cat("Iran cluster:", iran_cluster[1,23,1])Iran cluster: 3str(clusters)List of 9
$ cluster : int [1:223] 3 3 3 3 3 3 3 3 3 3 ...
$ centers : num [1:3, 1:21] 5.44e+11 1.52e+11 6.12e+09 5.85e+11 1.53e+11 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:3] "1" "2" "3"
.. ..$ : chr [1:21] "1996" "1997" "1998" "1999" ...
$ totss : num 9.08e+25
$ withinss : num [1:3] 6.21e+24 3.62e+24 3.11e+24
$ tot.withinss: num 1.29e+25
$ betweenss : num 7.79e+25
$ size : int [1:3] 2 12 209
$ iter : int 3
$ ifault : int 0
- attr(*, "class")= chr "kmeans"نتیجه فوق نشان می دهد ایران در دسته ی سوم قرار دارد.
۷. به وسیله تحلیل مولفه اصلی بعد داده رو به دو کاهش دهید سپس خوشه های به دست آمده در قسمت قبل را بر روی آن نمایش دهید. آیا عملکرد روش خوشه بندی شما مطلوب بوده است؟
ابتدا توسط pca بعد داده ها را به دو بعد کاهش می دهیم و سپس با kmean خوشه بندی می کنیم.
pca = prcomp(world_ec_cluster[,2:22], scale. = TRUE)
x = pca$x[,1:2]
world_ec_cluster <- world_ec_cluster %>% bind_cols(x %>% as.data.frame())
clusters = kmeans(world_ec_cluster[,24:25], 3)
world_ec_cluster$cluster_no = as.integer(clusters$cluster)
iran_cluster_pca = world_ec_cluster %>%
filter(country == 'Iran, Islamic Rep.')
cat("Iran cluster:", iran_cluster[1,23,1])Iran cluster: 3str(clusters)List of 9
$ cluster : int [1:223] 2 2 2 2 2 2 2 2 2 2 ...
$ centers : num [1:3, 1:2] 12.3806 -0.65044 53.17956 -0.8249 0.00992 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:3] "1" "2" "3"
.. ..$ : chr [1:2] "PC1" "PC2"
$ totss : num 4653
$ withinss : num [1:3] 352 292 0
$ tot.withinss: num 643
$ betweenss : num 4010
$ size : int [1:3] 7 215 1
$ iter : int 3
$ ifault : int 0
- attr(*, "class")= chr "kmeans"همانطور که مشاهده می کنیم هر دو دسته ی سوم را نشان می دهند. پس عملکرد خوشه بندی ما مناسب بوده است.
۸. با استفاده از داده روشی برای پیش بینی رشد اقتصادی ایران در سال آینده ارائه دهید.
از مدل glm استفاده می کنیم و سعی میکنیم سال ۲۰۱۶ را از روی سالهای پیشین بدست بیاوریم که نتایج نشان می دهد که به کمک داده های سال ۲۰۰۷ می توان پیشبینی کرد. اما همانطور که می بینیم مدل ما قوی نیست و خطای زیادی دارد.
iran_ec_m <- world_ec %>% filter(country == 'Iran, Islamic Rep.')
iran_ec_m[is.na(iran_ec_m)] = 0
glm_irn = glm(data = iran_ec_m[,3:23], formula = `2016`~.)
summary(glm_irn)
Call:
glm(formula = `2016` ~ ., data = iran_ec_m[, 3:23])
Deviance Residuals:
1 2 3 4 5 6 7 8
0.0001 0.1824 0.0002 0.0201 0.0001 0.0001 0.0001 2.6384
9 10 11 12 13 14 15 16
-0.8451 -3.1369 -1.0393 3.3487 -2.0293 2.5250 -3.1763 4.0549
17 18 19 20
-6.6281 4.0859 -0.0002 0.0002
Coefficients: (5 not defined because of singularities)
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2.1592 4.9213 -0.439 0.68351
`1996` 0.5198 0.7762 0.670 0.53971
`1997` 1.1219 1.3591 0.825 0.45550
`1998` -0.2441 0.9826 -0.248 0.81601
`1999` -0.8813 2.2597 -0.390 0.71641
`2000` -0.6201 0.6105 -1.016 0.36721
`2001` -0.3470 0.4529 -0.766 0.48627
`2002` -0.3611 1.5627 -0.231 0.82859
`2003` 1.9227 3.5487 0.542 0.61672
`2004` -1.0933 2.0732 -0.527 0.62584
`2005` -0.5921 0.9161 -0.646 0.55329
`2006` 1.2319 1.6720 0.737 0.50214
`2007` 0.4937 1.2457 0.396 0.71211
`2008` -1.8793 0.4548 -4.132 0.01447 *
`2009` NA NA NA NA
`2010` 0.1897 0.4410 0.430 0.68931
`2011` NA NA NA NA
`2012` NA NA NA NA
`2013` NA NA NA NA
`2014` 1.6837 0.3071 5.483 0.00539 **
`2015` NA NA NA NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for gaussian family taken to be 31.87369)
Null deviance: 2.5405e+24 on 19 degrees of freedom
Residual deviance: 1.2749e+02 on 4 degrees of freedom
AIC: 127.8
Number of Fisher Scoring iterations: 2glm_irn = glm(data = iran_ec_m[,3:23], formula = `2016`~`2007` + `2014`)
summary(glm_irn)
Call:
glm(formula = `2016` ~ `2007` + `2014`, data = iran_ec_m[, 3:23])
Deviance Residuals:
Min 1Q Median 3Q Max
-4.453e+10 4.140e+09 4.140e+09 4.140e+09 1.694e+10
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -4.140e+09 3.790e+09 -1.092 0.2900
`2007` 1.661e+00 7.209e-01 2.303 0.0341 *
`2014` -2.607e-01 6.048e-01 -0.431 0.6718
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for gaussian family taken to be 2.615323e+20)
Null deviance: 2.5405e+24 on 19 degrees of freedom
Residual deviance: 4.4460e+21 on 17 degrees of freedom
AIC: 1001.8
Number of Fisher Scoring iterations: 2glm_irn = glm(data = iran_ec_m[,3:23], formula = `2016`~`2007`)
summary(glm_irn)
Call:
glm(formula = `2016` ~ `2007`, data = iran_ec_m[, 3:23])
Deviance Residuals:
Min 1Q Median 3Q Max
-4.922e+10 4.057e+09 4.057e+09 4.057e+09 1.790e+10
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -4.057e+09 3.699e+09 -1.097 0.287
`2007` 1.350e+00 1.339e-02 100.777 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for gaussian family taken to be 2.497034e+20)
Null deviance: 2.5405e+24 on 19 degrees of freedom
Residual deviance: 4.4947e+21 on 18 degrees of freedom
AIC: 999.99
Number of Fisher Scoring iterations: 2lm_irn = lm(data = iran_ec_m[,3:23], formula = `2016`~.)
summary(lm_irn)
Call:
lm(formula = `2016` ~ ., data = iran_ec_m[, 3:23])
Residuals:
Min 1Q Median 3Q Max
-2010.82 0.00 43.97 233.67 266.80
Coefficients: (13 not defined because of singularities)
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2.683e+02 1.708e+02 -1.570 0.142
`1996` 1.257e+01 3.497e-02 359.510 <2e-16 ***
`1997` 4.314e+00 1.220e-03 3535.540 <2e-16 ***
`1998` 1.204e+00 4.457e-03 270.164 <2e-16 ***
`1999` -2.882e+01 7.615e-02 -378.413 <2e-16 ***
`2000` 7.695e+00 2.990e-02 257.326 <2e-16 ***
`2001` 5.642e+00 1.198e-02 470.847 <2e-16 ***
`2002` NA NA NA NA
`2003` NA NA NA NA
`2004` NA NA NA NA
`2005` 1.346e-03 3.103e-03 0.434 0.672
`2006` NA NA NA NA
`2007` NA NA NA NA
`2008` NA NA NA NA
`2009` NA NA NA NA
`2010` NA NA NA NA
`2011` NA NA NA NA
`2012` NA NA NA NA
`2013` NA NA NA NA
`2014` NA NA NA NA
`2015` NA NA NA NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 616.4 on 12 degrees of freedom
Multiple R-squared: 1, Adjusted R-squared: 1
F-statistic: 9.552e+17 on 7 and 12 DF, p-value: < 2.2e-16۹. سوالهای ۵ تا ۷ را ابتدا برای ۲۰ شاخص سلامت سپس بر حسب ۲۰ شاخص آموزشی تکرار کنید.
بهداشت
همانند بالا شاخص هایی که NA کمتری دارند را انتخاب می کنیم. در نهایت می بینیم ایران از نظر بهداشت نیز در دسته ی سوم قرار دارد. هم چنین در این بخش مدل خطی ما برای پیشبینی استفاده می کنیم. می بینیم خطای ما حول صفر است.
# health
# 5
world_hc <- wdi_data %>% filter(`Indicator Code`== 'SH.DYN.AIDS.FE.ZS') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]')) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'SH.DYN.AIDS.ZS') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]'))) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'SH.HIV.1524.MA.ZS') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]'))) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'SH.HIV.1524.FE.ZS') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]'))) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'SH.ANM.ALLW.ZS') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]'))) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'SH.PRG.ANEM') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]'))) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'SH.ANM.NPRG.ZS') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]'))) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'SH.ANM.CHLD.ZS') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]'))) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'SH.DTH.MORT') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]'))) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'SH.DTH.IMRT') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]'))) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'SH.VAC.TTNS.ZS') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]'))) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'SH.DYN.MORT') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]'))) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'SH.DYN.NMRT') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]'))) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'SH.STA.MMRT') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]'))) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'SH.MMR.RISK.ZS') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]'))) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'SH.XPD.GHED.PP.CD') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]'))) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'SH.XPD.GHED.PC.CD') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]'))) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'SH.XPD.CHEX.PC.CD') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]'))) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'SH.HIV.INCD') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]'))) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'SH.HIV.0014') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]')))
world_hc_m <- world_hc %>%
melt(id.vars=c('country', 'index')) %>%
filter(!is.na(value)) %>%
group_by(index)
plots <- world_hc_m %>%
do(
plots = ggplot(.,aes(x = variable, y = value, color = index)) + geom_boxplot() +
xlab("Year") + ggtitle("Worldwide Health vs. Iran") +
theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
geom_line(data = subset(.,country == 'Iran, Islamic Rep.'),aes(x = variable, y = value, group = 1, color = 'Iran, Islamic Rep.'))
)
plots$plots[[1]]
[[2]]
[[3]]
[[4]]
[[5]]
[[6]]
[[7]]
[[8]]
[[9]]
[[10]]
[[11]]
[[12]]
[[13]]
[[14]]
[[15]]
[[16]]
[[17]]
[[18]]
[[19]]
[[20]]
# 6
# cleaning data
world_hc_cluster <- world_hc
world_hc_cluster[is.na(world_hc)] = 0
patt = "Africa|America|Asia|dividend|Euro|income|members|Middle|only|Sub|total|World"
world_hc_cluster <- world_hc_cluster %>% group_by(country) %>%
filter(!grepl(patt,country)) %>%
summarise_at(c(2:22), mean)
set.seed(134)
clusters = kmeans(world_hc_cluster[,2:22], 3)
world_hc_cluster$cluster_no = as.integer(clusters$cluster)
iran_cluster = world_hc_cluster %>%
filter(country == 'Iran, Islamic Rep.')
cat("Iran cluster:", iran_cluster[1,23,1])Iran cluster: 3str(clusters)List of 9
$ cluster : int [1:223] 3 3 3 3 3 3 3 3 3 3 ...
$ centers : num [1:3, 1:21] 41381 215984 1876 40758 214704 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:3] "1" "2" "3"
.. ..$ : chr [1:21] "1996" "1997" "1998" "1999" ...
$ totss : num 3.76e+12
$ withinss : num [1:3] 1.00e+11 2.18e+11 3.41e+10
$ tot.withinss: num 3.52e+11
$ betweenss : num 3.41e+12
$ size : int [1:3] 11 5 207
$ iter : int 3
$ ifault : int 0
- attr(*, "class")= chr "kmeans"# 7
pca = prcomp(world_hc_cluster[,2:22], scale. = TRUE)
x = pca$x[,1:2]
world_hc_cluster <- world_hc_cluster %>% bind_cols(x %>% as.data.frame())
clusters = kmeans(world_hc_cluster[,24:25], 3)
world_hc_cluster$cluster_no = as.integer(clusters$cluster)
iran_cluster_pca = world_hc_cluster %>%
filter(country == 'Iran, Islamic Rep.')
cat("Iran cluster:", iran_cluster[1,23,1])Iran cluster: 3str(clusters)List of 9
$ cluster : int [1:223] 2 2 2 2 2 2 2 2 2 2 ...
$ centers : num [1:3, 1:2] -4.36247 0.90214 -27.75121 0.05342 -0.00315 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:3] "1" "2" "3"
.. ..$ : chr [1:2] "PC1" "PC2"
$ totss : num 4661
$ withinss : num [1:3] 132.4 42.9 257.4
$ tot.withinss: num 433
$ betweenss : num 4228
$ size : int [1:3] 11 207 5
$ iter : int 3
$ ifault : int 0
- attr(*, "class")= chr "kmeans"# 8
iran_hc_m <- world_hc %>% filter(country == 'Iran, Islamic Rep.')
iran_hc_m[is.na(iran_hc_m)] = 0
lm_irn = lm(data = iran_hc_m[,3:23], formula = `2016`~.)
summary(lm_irn)
Call:
lm(formula = `2016` ~ ., data = iran_hc_m[, 3:23])
Residuals:
1 2 3 4 5 6
-1.161e-19 8.599e-10 -4.793e-10 -3.806e-10 2.435e-19 -2.661e-20
7 8 9 10 11 12
-1.428e-19 3.862e-20 2.169e-22 -3.504e-23 -5.071e-22 7.068e-20
13 14 15 16 17 18
-4.862e-19 -3.858e-20 -5.987e-18 -2.820e-22 1.192e-21 -3.030e-22
19 20
-2.412e-23 -1.656e-22
Coefficients: (3 not defined because of singularities)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.081e-02 4.460e-10 4.666e+07 4.59e-16 ***
`1996` 1.606e+01 7.053e-08 2.277e+08 < 2e-16 ***
`1997` -6.925e+01 3.826e-07 -1.810e+08 < 2e-16 ***
`1998` 1.650e+02 9.457e-07 1.745e+08 < 2e-16 ***
`1999` -1.290e+02 7.238e-07 -1.782e+08 < 2e-16 ***
`2000` -4.657e+01 1.808e-07 -2.576e+08 < 2e-16 ***
`2001` -2.141e+01 2.585e-07 -8.283e+07 < 2e-16 ***
`2002` 8.476e+01 5.288e-07 1.603e+08 < 2e-16 ***
`2003` 1.032e+02 5.522e-07 1.869e+08 < 2e-16 ***
`2004` -1.125e+02 6.232e-07 -1.806e+08 < 2e-16 ***
`2005` 3.327e+01 1.634e-07 2.036e+08 < 2e-16 ***
`2006` -6.425e+01 2.408e-07 -2.668e+08 < 2e-16 ***
`2007` 3.470e+01 2.807e-08 1.236e+09 < 2e-16 ***
`2008` 5.312e+01 3.681e-07 1.443e+08 < 2e-16 ***
`2009` -6.090e+01 2.908e-07 -2.094e+08 < 2e-16 ***
`2010` 2.240e+01 8.453e-08 2.650e+08 < 2e-16 ***
`2011` -1.838e+01 1.089e-08 -1.688e+09 < 2e-16 ***
`2012` 1.058e+01 4.936e-08 2.143e+08 < 2e-16 ***
`2013` NA NA NA NA
`2014` NA NA NA NA
`2015` NA NA NA NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 7.463e-10 on 2 degrees of freedom
Multiple R-squared: 1, Adjusted R-squared: 1
F-statistic: 6.783e+25 on 17 and 2 DF, p-value: < 2.2e-16
آموزش
همانند بالا شاخص هایی که NA کمتری دارند را انتخاب می کنیم. در نهایت می بینیم ایران از نظر آموزش نیز در دسته ی سوم قرار دارد. هم چنین در این بخش مدل خطی ما برای پیشبینی استفاده می کنیم. می بینیم خطای ما حول صفر است.
# education
# 5
world_ed <- wdi_data %>% filter(`Indicator Code`== 'SE.ENR.PRIM.FM.ZS') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]')) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'SE.PRM.DURS') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]'))) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'SE.XPD.TERT.ZS') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]'))) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'SE.XPD.PRIM.Z') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]'))) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'SE.XPD.SECO.ZS') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]'))) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'SE.PRM.ENRR') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]'))) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'SE.PRM.TCHR') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]'))) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'SE.PRM.ENRL') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]'))) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'SE.PRE.ENRL.TC.ZS') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]'))) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'SE.PRE.ENRR') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]'))) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'SE.SEC.PROG.ZS') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]'))) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'SE.SEC.ENRL.GC') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]'))) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'SE.PRM.PRSL.ZS') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]'))) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'SE.SEC.AGES') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]'))) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'SE.PRE.ENRL.TC.ZS') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]'))) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'SE.PRM.TENR.FE') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]'))) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'SE.SEC.ENRL') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]'))) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'SE.PRM.NENR') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]'))) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'SE.SEC.DURS') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]'))) %>%
rbind(wdi_data %>% filter(`Indicator Code`== 'SE.PRM.AGES') %>%
select(country = `Country Name`, index = `Indicator Name`, matches('^199[6-9]|^20[0-1][0-9]')))
world_ed_m <- world_ed %>%
melt(id.vars=c('country', 'index')) %>%
filter(!is.na(value)) %>%
group_by(index)
plots <- world_ed_m %>%
do(
plots = ggplot(.,aes(x = variable, y = value, color = index)) + geom_boxplot() +
xlab("Year") + ggtitle("Worldwide Education vs. Iran") +
theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
geom_line(data = subset(.,country == 'Iran, Islamic Rep.'),aes(x = variable, y = value, group = 1, color = 'Iran, Islamic Rep.'))
)
plots$plots[[1]]
[[2]]
[[3]]
[[4]]
[[5]]
[[6]]
[[7]]
[[8]]
[[9]]
[[10]]
[[11]]
[[12]]
[[13]]
[[14]]
[[15]]
[[16]]
[[17]]
[[18]]
# 6
# cleaning data
world_ed_cluster <- world_ed
world_ed_cluster[is.na(world_ed)] = 0
patt = "Africa|America|Asia|dividend|Euro|income|members|Middle|only|Sub|total|World"
world_ed_cluster <- world_ed_cluster %>% group_by(country) %>%
filter(!grepl(patt,country)) %>%
summarise_at(c(2:22), mean)
set.seed(134)
clusters = kmeans(world_ed_cluster[,2:22], 3)
world_ed_cluster$cluster_no = as.integer(clusters$cluster)
iran_cluster = world_ed_cluster %>%
filter(country == 'Iran, Islamic Rep.')
cat("Iran cluster:", iran_cluster[1,23,1])Iran cluster: 1str(clusters)List of 9
$ cluster : int [1:223] 2 2 2 2 2 2 2 2 2 2 ...
$ centers : num [1:3, 1:21] 1094749 77381 13089401 923136 66614 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:3] "1" "2" "3"
.. ..$ : chr [1:21] "1996" "1997" "1998" "1999" ...
$ totss : num 8.97e+15
$ withinss : num [1:3] 3.37e+14 1.37e+14 2.61e+14
$ tot.withinss: num 7.35e+14
$ betweenss : num 8.23e+15
$ size : int [1:3] 15 206 2
$ iter : int 3
$ ifault : int 0
- attr(*, "class")= chr "kmeans"# 7
pca = prcomp(world_ed_cluster[,2:22], scale. = TRUE)
x = pca$x[,1:2]
world_ed_cluster <- world_ed_cluster %>% bind_cols(x %>% as.data.frame())
clusters = kmeans(world_ed_cluster[,24:25], 3)
world_ed_cluster$cluster_no = as.integer(clusters$cluster)
iran_cluster_pca = world_ed_cluster %>%
filter(country == 'Iran, Islamic Rep.')
cat("Iran cluster:", iran_cluster[1,23,1])Iran cluster: 1str(clusters)List of 9
$ cluster : int [1:223] 2 2 2 2 2 2 2 2 2 2 ...
$ centers : num [1:3, 1:2] 42.8128 -0.7341 4.3726 -2.8555 -0.0809 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:3] "1" "2" "3"
.. ..$ : chr [1:2] "PC1" "PC2"
$ totss : num 4547
$ withinss : num [1:3] 206.2 70.1 155.7
$ tot.withinss: num 432
$ betweenss : num 4115
$ size : int [1:3] 2 206 15
$ iter : int 3
$ ifault : int 0
- attr(*, "class")= chr "kmeans"# 8
iran_ed_m <- world_ed %>% filter(country == 'Iran, Islamic Rep.')
iran_ed_m[is.na(iran_ed_m)] = 0
lm_irn = lm(data = iran_ed_m[,3:23], formula = `2016`~.)
summary(lm_irn)
Call:
lm(formula = `2016` ~ ., data = iran_ed_m[, 3:23])
Residuals:
1 2 3 4 5 6
-4.361e-01 4.475e-01 -1.263e-04 5.485e-06 -1.001e-04 9.614e-09
7 8 9 10 11 12
-1.343e-08 -8.379e-05 -4.014e-04 -1.630e-05 2.662e-11 1.349e-06
13 14 15 16 17 18
-3.400e-01 -8.379e-05 3.342e-05 -1.661e-10 -4.751e-05 7.039e-02
19
2.590e-01
Coefficients: (5 not defined because of singularities)
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.04650 0.42536 -2.460 0.090843 .
`1996` 0.51980 0.18541 2.803 0.067658 .
`1997` 0.22063 0.01284 17.185 0.000429 ***
`1998` 1.44937 0.14168 10.230 0.001991 **
`1999` 9.12301 0.47831 19.074 0.000315 ***
`2000` -12.98219 0.62537 -20.759 0.000244 ***
`2001` -56.82145 2.92022 -19.458 0.000297 ***
`2002` 19.89924 0.97942 20.317 0.000261 ***
`2003` 130.84586 6.72701 19.451 0.000297 ***
`2004` -99.89070 5.02793 -19.867 0.000279 ***
`2005` 23.11709 1.08352 21.335 0.000225 ***
`2006` -3.19842 0.31194 -10.253 0.001978 **
`2007` -62.44794 3.16876 -19.707 0.000285 ***
`2008` 49.74231 2.55809 19.445 0.000297 ***
`2009` NA NA NA NA
`2010` 0.61198 0.03251 18.827 0.000327 ***
`2011` NA NA NA NA
`2012` NA NA NA NA
`2013` 0.94268 0.27009 3.490 0.039762 *
`2014` NA NA NA NA
`2015` NA NA NA NA
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.439 on 3 degrees of freedom
Multiple R-squared: 0.9997, Adjusted R-squared: 0.9983
F-statistic: 687.4 on 15 and 3 DF, p-value: 8.032e-05۱۰. کشورهای دنیا را بر حسب ۶۰ شاخص اقتصادی، سلامت و آموزش با روش سلسله مراتبی خوشه بندی کرده و دندروگرام آن را رسم نمایید. اگر داده ها بر سه دسته تقسیم شوند ایران در کدام دسته می گنجد؟
۱۱. سه یافته جالب از داده ها استخراج کنید.